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We study the Hubbard model on a bilayer with repulsive on-site interactions, U, in which fermions 
undergo both intra-plane (f) and inter-plane (t z ) hopping. This situation is what one would expect 
in high-temperature superconductors such as YBCO, with two adjacent Cu02 planes. Magnetic 
and pairing properties of the system are investigated through Quantum Monte Carlo simulations 
for both half- and quarter-filled bands. We find that in all cases inter-planar pairing with d x i_ z i 
symmetry is dominant over planar pairing with d x 2_ y 2 symmetry, and that for t z large enough pair 
formation is possible through antiferromagnetic correlations. However, another mechanism is needed 
to make these pairs condense into a superconducting state at lower temperatures. We identify the 
temperature for pair formation with the spin gap crossover temperature. 
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I. INTRODUCTION 

The normal state of the high temperature supercon- 
ductors unveils interesting features, especially in the 
YBCO compounds. Firstly, the in-plane resistivity of 
optimally doped (i.e., hole fraction x = x m , correspond- 
ing to a maximum T c ) samples of YB^CuaOy-j, with 
y = 0.10, displays linear behavior in a wide range of 
temperatures, from just above T c ~ 90 K to nearly 1000 
K; 1 ' 2 this is to be contrasted with the T 2 behavior at low 
temperatures, that would be expected for a usual Fermi 
liquid. For underdoped samples (i.e., x < x m ), the linear 
temperature behavior of the resistivity crosses over to T a , 
with a ~ 2.5, below a characteristic temperature T x (a;). 3 
Overdoped samples, on the other hand, behave like usual 
metals, in the sense that the resistivity crosses over from 
a linear-T dependence (at high temperatures) to Fermi 
liquid behavior at lower temperatures. 4 ' 5 Secondly, the 
magnetic response - probed by nuclear relaxation rate 
(1/TT) on the Cu(2) sites 6 ' 7 and by NMR Knight shift 
at the 89 Y sites 8 - also shows different behavior in the 
underdoped and overdoped regimes. For overdoped sam- 
ples, the relaxation rate, which is proportional to the 
staggered susceptibility (x(k)), increases monotonically 
as the temperature is lowered, in a Curie-like fashion; 6 ' 7 
this is indicative of gapless spin excitations from a state 
with strong antiferromagnetic correlations. At T c one 
expects the simultaneous formation and condensation of 
pairs, leading, respectively, to spin and superconduct- 
ing gaps. The situation is similar for the NMR Knight 
shift, which is proportional to the uniform susceptibility 
(x(0)): one observes a Pauli-like behavior, in the sense 
that it is roughly temperature independent. 3 For under- 
doped samples both xi^) an d x(0) show a marked devi- 
ation from their respective Curie and Pauli behaviors as 
the temperature is lowered; i.e., they start to decrease 
as T decreases below certain temperature scales T^(x) 
and To (a;). While one cannot be precise to the point of 
identifying all temperature scales, the general expecta- 



tion is that T x ~ T T ~ To. Thus, one may interpret 
the behavior in the underdoped region as due to the for- 
mation of pairs (without condensation) at a temperature 
T x , which is accompanied by the opening of a spin gap. 9 
As the temperature is lowered further, a superconducting 
gap opens at T c , giving rise to a superconducting insta- 
bility. Further, Ito et al. 3 have stressed that the temper- 
ature dependence of their in-plane resistivity data, p(T), 
is such that x(0) ~ p(T)/T and x(tt) ~ p(T)/T 2 . Similar 
conclusions may be drawn from the analysis of transport 
properties of YBa2Cu40s, corresponding to the under- 
doped regime. 10 The interpretation of these data is still a 
matter of debate. On the one hand, the above mentioned 
transport (charge) properties of the underdoped materi- 
als are strongly influenced by the opening of a spin gap, 
suggesting inseparability of spin and charge behaviors. 11 
On the other hand, one could also assume spin and charge 
separation in the sense that the spin gap opening would 
indicate formation of pairs, which would Bose condense 
(superconduct) at a lower temperature; 9 ' 12 the overdoped 
regime would then be described as an ordinary Fermi liq- 
uid, with both spin and superconducting gaps opening at 
the same critical temperature, T c . 

In order to gain insight into these issues, one works 
with simplified models which may highlight the main 
physical mechanisms at work. For instance, one possi- 
bility is to treat the planar spin excitations as those of a 
nearly antiferromagnetic Fermi liquid (NAFL). 13 ' 14 This 
approach has been used to discuss several properties of 
overdoped compounds, 13-16 though it is not clear at the 
moment how to incorporate the spin gap in the treat- 
ment of the underdoped regime. An alternative explana- 
tion for the experimental data is based on the presence 
of a double Cu02 layer in YBCO compounds, as oposed 
to the single layered structure of the La2- a: Sr ;l ;Cu04 
materials. 17 While the spin gap behavior in underdoped 
La compounds was attributed 17 to a spin-density wave 
(SDW) instability, it was suggested that antiferromag- 
netic correlations between fermions in adjacent planes 
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in bilayer materials lead to singlet pairing with each 
member of the pair lying on each plane; this, in turn, 
would be responsible for the non-Fermi liquid behavior. 
These ideas 17-19 have been tested on models (e.g., the t- J 
model) with antiferromagnetic Heisenberg-like coupling 
between spins on different layers, 12 ' 17 ' 20 ' 21 and they in- 
deed lead to pairing between fermions in different planes. 
Unfortunately, interplane couplings of the order of (and, 
in some cases, even larger than) intraplane couplings are 
generally needed to achieve pairing, which seems rather 
unrealistic. The Hubbard model on a bilayer, with on- 
site repulsion and both intra- and inter-plane hopping, 
can be considered as a weaker- coupling version of the 
models mentioned above, and should therefore provide a 
more realistic description of the actual physical situation. 
The Hamiltonian is 

n = - Uj (c] a c jtT + H.c.) 

(iJ) 

where the sums run over sites of two square lattice lay- 
ers, denotes nearest neighbor sites, and the hopping 
integral is given by 

_ | t if i and j are within the same plane , , 
ZJ \ t z if i and j lie in different planes . ^ ' 

c lo- ( c i<r) creates (annihilates) a fermion at site i with spin 
<t, H.c. stands for Hermitian conjugate, U > is the on- 
site repulsion, and /i is the chemical potential controlling 
the band-filling, (n). This model has been studied pre- 
viously by Bulut et al. 22 , who found evidence for node- 
less singlet pairing, from random-phase-approximation- 
Eliashberg calculations based on the exchange of spin 
fluctuations, for (n) = 1 and 0.85; their results for half- 
filling were compared with those from Monte Carlo sim- 
ulations on a 2x(4x4) lattice. In view of the impor- 
tant role interlayer pairing may play as a mechanism for 
high-temperature superconductivity, independent checks 
should be made in order to assure the effect is indepen- 
dent of model details and of the approximations used. 
With this in mind, here we report the results of extensive 
Monte Carlo simulations, in which both magnetic and 
pairing properties of the model are examined for larger 
lattices [up to 2x(8x8)] and for the cases of half- and 
quarter-filled bands. 

The layout of the paper is as follows. In Sec. II we 
mention the difficulties posed by the "minus sign prob- 
lem" and present the results for magnetic susceptibili- 
ties and correlation functions. Similarly, Sec. Ill deals 
with the superconducting susceptibilities and correlation 
functions. Finally, Sec. IV summarizes our findings and 
presents the conclusions. 



II. MAGNETIC PROPERTIES 

In a grand-canonical quantum Monte Carlo 
simulation 23-27 the imaginary time is discretized through 
the introduction of M "time" slices, with M <x l/T. The 
"Boltzmann weight" is given by a product of fermion 
determinants, and is only positive definite for the repul- 
sive model at half-filling and for the attractive model 
at any filling; 24 ' 27 ' 28 otherwise, some configurations will 
give negative contribution. This "minus sign problem" 
is circumvented by redefining the averages in such way 
that an average sign appears in the denominator. 29 The 
average sign, (e)', behaves generically as follows: 30 for a 
fixed temperature, it drops dramatically from 1 very near 
half-filling, reaching a minimum for some band filling, 
and eventually increasing to (e) ~ 1 for some occupation 
near (n) = 1/2; the value of (n) for which this happens is 
not very sensitive to details such as the values of U, of t z , 
and so on. 30 As one dopes further away, the average sign 
drops again. Moreover, as the temperature decreases, 
the dip in (e) becomes more pronounced, reaching a very 
small value; this introduces uncontrollable noise in the 
calculation of average values. In view of this, we restrict 
our discussion to half- and quarter- filled bands (i.e., to 
(n) = 1 and f/2, respectively). 

The clusters used here have N s = 2 x (L x L) sites, 
with periodic boundary conditions (PBC); that is, each 
site is connected with its six nearest neighbors through 
a hopping term. Due to the PBC along a direction of 
finite size L z = 2, the effective hopping along the z direc- 
tion becomes tj_ = 2t z . The simulations were performed 
on Sun and IBM RISC 600/525 workstations, and on a 
CRAY Y-MP/2E; a single datum point involves between 
500 and 2000 MC sweeps over all time slices and we took 
At = (3/M= 0.125. 

Magnetic properties are examined through both the 
magnetic susceptibility, 

= W s E e!q ' (r, ~ rj) f dT K0>;(°)> > ( 3 ) 

with 

m;(T) = e T *[n it -n i Je- T * , (4) 
and the magnetic structure factor, 

s(q) = ^-E e!q ' (r, ~ rj) <™i™;>. ( 5 ) 

where m; = rrii(0). 

Figure f shows the uniform (q = 0) and staggered 
(q = 77) susceptibilities for the half-filled band case, with 
tj_ = 0.7 and [7 = 4 (from now on, energies will be ex- 
pressed in units of t, the nearest-neighbor planar hopping 
integral, and temperatures in units of t/kB, ks being the 
Boltzmann constant); the system sizes are 2 X (L X L). 
with L = 4, 6 and 8. The uniform susceptibility (see Fig. 
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1(a)) decreases as the temperature decreases, in a non- 
Pauli fashion; this should be contrasted with the behavior 
of the staggered susceptibility (see Fig. 1(b); notice the 
difference in scales), which shows a Curie- like behavior 
and scaling with lattice size, L. By analogy with the 
single-layer Hubbard model, 24 ' 31 one associates an insu- 
lating antiferromagnetic state also for the bilayer at half 
filling. 



maxima and roughly the same heights, though one may 
argue that as tj_ increases, the maximum in the q z = 
sector narrows slightly, while that in the q z = 7r sector 
flattens out also slightly. The important point is that 
antiferromagnetic correlations are always dominant, at 
least in one (planar) direction, signalled by <S*(q) being 
roughly the same, as long as q x = 77. 
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FIG. 1. Uniform (a) and staggered (b) magnetic suscepti- 
bilities as functions of temperature, for the repulsive Hubbard 
bilayer at half-filling, with U = 4 and t± = 0.7. System sizes 
are 2 X (L X L), with L = 4 (triangles), L = 6 (squares), and 
L = 8 (circles). The lines are guides to the eye. 

As mentioned before, the minus-sign problem prevents 
us from going to very low temperatures away from half 
filling; the optimum choice for the bilayers i.e., the one 
that allows us to reach the lowest temperatures with an 
acceptable (e) > 0.6, is (n) = 0.5. In order to discuss the 
magnetic correlations for this filling, in Fig. 2 we display 
the structure factor, Eq. (5), at fixed temperature; since 
the data are symmetric under the exchange (q x ,q y ) — ^ 
(q y ,q x ), the Brillouin zone is only shown partially. In 
each of the q z sectors the behavior is similar for tj_ = 0.2 
[Fig. 2(a)] and for t± = 0.7 [Fig. 2(b)], with very broad 
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FIG. 2. Wavevector dependence of the structure factor for 
a fixed inverse temperature, /3 = 8, for the repulsive Hubbard 
bilayer at quarter-filling, with U = 4 and t± = 0.2 (a) and 
fx = 0.7 (b). System sizes are 2 X (L X L), with L = 6 
(squares), and L = 8 (circles). The error bars are smaller 
than the data points, and the lines are guides to the eye. 

It is also instructive to discuss the q-dependence of 
the susceptibility. For tj_ = 0.2 (see Fig. 3(a)), two 
peaks with approximately the same height appear in both 
q z = and q z = tt sectors, indicating that the magnetic 
response is the same whether the field is uniform or stag- 
gered along the z direction (perpendicular to the two lay- 
ers). As one increases the interlayer hopping to tj_ = 0.7 
(see Fig. 3(b)), only one peak is found, corresponding to 
a field staggered in the z-direction. These peaks come 
about as a result of correlations adding up coherently; 
we shall return to this point below. 
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FIG. 3. Wavevector dependence of the susceptibility for a 
fixed inverse temperature, /3 = 8, for the repulsive Hubbard 
bilayer at quarter-filling, with U = 4 and t± = 0.2 (a) and 
fx = 0.7 (b). System sizes are 2 X (L X L), with i = 6 
(squares), and L = 8 (circles). The lines are guides to the 
eye. 

Figure 4(a) shows the temperature dependence of the 
uniform susceptibility at quarter-filled band for tj_ = 0.2, 
and for L = 6 and 8. From the data for L = 6 alone, 
one might be tempted to infer a spin gap behavior, sig- 
nalled by the downturn of x(0), together with the absence 
of long range magnetic order. This behavior, however, 
does not seem to persist for the L = 8 system; on the 
contrary, the weak temperature dependence of x(0) sug- 
gests a Pauli-like — and therefore metallic — behavior 
for the whole temperature range examined. The data for 
tj_ = 0.7, shown in Fig. 4(b), can be interpreted differ- 
ently. As the temperature decreases, a downturn in x(0) 
is followed by an increase at lower temperatures; unlike 
what was observed in Fig. 4(a), this feature is indepen- 
dent of system size. While this is not the prototypical 
spin gap behavior, the different temperature dependences 
found for the two values of t± can hardly be regarded 
as fortuitous. This, together with the above analysis of 
x(q), suggests that the nature of spin excitations changes 



as one increases tj_. In the following section we discuss 
the possible bearings of the magnetic properties on pair- 
ing. 
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FIG. 4. Uniform susceptibility as a function of tempera- 
ture, for the repulsive Hubbard bilayer at quarter-filling, with 
U = 4 and t± = 0.2 (a) and t± = 0.7 (b). System sizes are 
2 X (L X L), with L = 6 (squares), and L = 8 (circles). The 
lines are guides to the eye. 



III. SUPERCONDUCTING PROPERTIES 

Superconducting properties are probed by the uniform 
(q = 0) zero-frequency pair susceptibilities, 

P x = fdr <A A (r)A+(0)) , (6) 
Jo 

and by equal-time uniform (q = 0) correlation functions, 
C x = (A{A X + A X A{) , (7) 
where the pair-field operator is given by 

= ^Ew k )4A, (8) 
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with /A(k) defining its symmetry. It is by now well 
established 32-35 that the dominant pairing correlations 
in the planar Hubbard model have d x i_ y i -symmetry, for 
which 



fd^_, (k) = cosk x - cos k y 



(9) 



Our strategy here is to compare this dominant planar 
pairing with inter-planar correlations of different sym- 
metries, such as 



fd z (k) = cos k z , 
fd x 2_,2 (k) = cos k x - cos k z , 



f s , (k) = cos k x + cos k z 



fd x ,(k) = s'mk x s'mk z 



(10) 



(11) 



(12) 



(13) 



as well as with other combinations of longer range. 

As it turned out, the dominant susceptibility corre- 
sponding to interplanar pairing has d x 2_ z 2 symmetry; 
that is, it is larger than any other in all situations (i.e., 
different temperatures, interlayer hoppings, and band fill- 
ings) examined, especially the one corresponding to d z 
(the so-called nodeless d-wave 22 ). We have also checked 
that the largest planar pairing susceptibility has d x i_ y i 
symmetry even in the case of a bilayer. In what follows, 
we therefore concentrate our discussion in terms of the 
dominant ones. 
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FIG. 5. Temperature dependence of the pairing suscepti- 
bilities at quarter-filling, for t± = 0.2, for different system 
sizes: 2 X (L X L), with (a) L = 4, (b) L = 6, and (c) L = 8. 
Data for the interacting case with U = 4 [non-interacting], are 
represented by squares [dotted] for d x 2_ z 2 symmetry, and by 
triangles [dashed] for d x 2_ y 2 symmetry; solid lines are guides 
to the eye. 
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FIG. 6. Same as Fig. 5, but for t± = 0.7 

Figure 5 shows the temperature dependence of Pj 2 2 
and Pd 2 2 5 l° r l ne interacting (U = 4) and free (U = 0) 
cases, for different system sizes, at quarter filling, and 
with tj_ = 0.2. In all cases, the susceptibilities are sup- 
pressed by the presence of a repulsive on-site interac- 
tion; this means that the fermions are less likely to pair 
in the presence of the on-site repulsion, similarly to the 
single-layer Hubbard model. 32 Nevertheless, it is impor- 
tant to note that interplanar pairing always dominates 
the planar one. As one increases the interlayer hopping, 
an interesting feature emerges. In Fig. 6, for tj_ = 0.7, 
the susceptibilities in the interacting case grow faster at 
lower temperatures than those corresponding to the free 
case, for the largest system examined. One therefore ex- 
pects that there is a temperature, T* ~ 0.1, below which 
the interacting susceptibility is larger than the free one, 
indicating a tendency towards pairing arising from repul- 
sion. In view of the analysis of Sec. II, one concludes that 
pairing is favored when "time" -correlations along the z 
direction are predominantly of antiferromagnetic nature; 
that is, when the degeneracy between uniform and stag- 
gered correlations is broken. 

In order to verify whether the system could be close 
to an actual phase transition, we have also examined the 
dependence of the pairing correlation function, Eq. (7), 
with the system size. If an infinite two-dimensional sys- 
tem (or bilayer) undergoes a superconducting transition, 
it belongs to the Kosterlitz-Thouless 36 XK-model uni- 
versality class. Accordingly, pairing correlations should 
become critical at a critical temperature T c , and decay 
algebrically, 



G(r) ~ r 



(14) 



with r) = 0.25 for T —¥ T c + . From finite-size scaling (FSS) 
theory, 37 one then infers that for a system of linear size 
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L, its associated uniform Fourier transform becomes 



/' 

Jo 



d 2 r r 



L 2- V 



(15) 



Above criticality, the apropriate scaling variable 37 is L/£, 
where £ ~ exp(^4/ y/T — T c ) , with A being of order unity, 
is the correlation length for the infinite system. There- 
fore, one can assume the following FSS ansatz: 

C x (T,L) = L^F(L/0, (16) 

where F(z) is a scaling function such that F(z) — ^ 
constant when L <C £, recovering Eq. (15). At T c , £ = oc. 
so that L V ~ 2 C\(T C , L) is a constant independent of lat- 
tice size. By plotting L n ~ 2 C\(T, L) as a function of T 
for systems of different sizes, a phase transition would 
be signalled by a crossing of curves corresponding to sys- 
tems of different size. 38 Figure 7 shows the size-scaled 
pairing correlation function with d x 2_ z 2 symmetry, as a 
function of the inverse temperature. For all sizes studied 
L~ 7 / 4 C stabilizes to a constant value at large /3, with- 
out any indication of crossing. The possibility of a phase 
transition is therefore ruled out, and the behavior of the 
pairing susceptibility should be attributed to pair for- 
mation. These preformed pairs should then condense — 
through a mechanism absent in the present model — at 
a lower temperature into a superconducting state. 
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FIG. 7. Inverse-temperature dependence of the size-scaled 
uniform i x i_ z i— pairing correlation function at quarter-filling, 
for t± = 0.7 and different system sizes: 2x(Ixi), with L = 4 
(triangles), L = 6 (squares), and L = 8 (circles). Solid lines 
are guides to the eye. 



IV. CONCLUSIONS 

We have examined the interplay between magnetism 
and pairing in the Hubbard model on a bilayer. We have 



found that, even though the magnetic structure factor 
is not sharply peaked at any particular wavevector, some 
correlations add up coherently, giving rise to peaks in the 
q-dependent susceptibility. For smaller interlayer hop- 
ping, the system's response is similar whether the field 
is staggered or uniform along the direction perpendicular 
(z) to the bilayer. For larger interlayer hopping this de- 
generacy is broken, and the system is only distinctively 
responsive to a z-staggered field. 

With respect to pairing, we have found that inter- 
layer correlations with d x 2_ z 2 symmetry are the domi- 
nant ones, in the sense that their associated susceptibilty 
is larger than any other, including planar ones. Thus, 
the suggestion by Milllis and Monien 17 has been con- 
firmed away from the strong-coupling regime. We have 
also found evidence that larger interlayer hopping (i.e., 
tj_ = 0.7) may lead to pairing at low temperatures, and 
associate this to the dominant antiferromagnetic ( "time"- 
) correlations between the two layers. On the other hand, 
a finite-size scaling analysis of the pairing correlations in- 
dicates the absence of a phase transition into a state with 
long range (or quasi-long range) order. Upon completion 
of this work, we have become aware of two recent studies. 
Hetzel et al. 39 discussed Hubbard bilayers at zero temper- 
ature, and have also found no evidence for off-diagonal 
long range order in the model. Scalettar et al. 40 have 
considered the half-filled case at finite temperatures, in- 
cluding different chemical potentials on each plane, and 
have also ruled out off-diagonal long-range order. 

The picture that emerges is that pair formation in Hub- 
bard bilayers is possible through antiferromagnetic cor- 
relations, though another mechanism is needed to make 
these pairs condense into a superconducting state at 
lower temperatures. In this respect, the situation is sim- 
ilar to what happens in the underdoped regime of YBCO 
compounds, as mentioned in the Introduction. The tem- 
perature T* , at which the interacting susceptibility be- 
comes larger than the non-interacting one should be as- 
sociated with T x , the spin-gap crossover temperature. 
Indeed, if one takes T* = 0.05t and typical values for 
the Hubbard model parameters (t ~ 0.5 eV, U ~ At ~ 
2 eV), we obtain the estimate T* ~ 250 K, which is 
of the correct order of magnitude for T x in YBCO. Fi- 
nally, we should mention that this picture has appealing 
similarities with the results obtained from the attractive 
Hubbard model; see e.g., Ref. 41. 
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